#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include "common.h"
#include <list>

using namespace std;
#define MAX 6   // Maximum number of particles possible in a box


// Inits a data structure 'boxes' equivalent to the following
// list<int> boxes[n][n]
// Provided the restriction that the maximum elements in a list can only be
// MAX, and no more. Internally, it is a 3D array of the form 
// int boxes[n][n][MAX+1], where boxes[x][y][0] indicates the current
// number of elements in the list boxes[x][y]
inline int* init_ds(int n) {
    return (int *)calloc(n*n*(MAX+1), sizeof(int));
}

// This cleans up the boxes (equivalent to calling clear on all lists)
inline void clean_ds(int *boxes, int n) {
    memset(boxes, 0, sizeof(int)*n*n*(MAX+1));
}

// This frees up the memory
inline void dest_ds(int *boxes) {
    free(boxes);
}

// find boxes[x][y][z]
inline int get_val(int *boxes, int n, int x, int y, int z) {
    return boxes[((x*n)+y)*(MAX+1) + z];
}

// set boxes[x][y][z] = n
inline void set_val(int *boxes, int n, int x, int y, int z, int val) {
    boxes[((x*n)+y)*(MAX+1) + z] = val;
}   

// Equivalent to inserting a value in the list boxes[x][y]
inline void insert_ds(int *boxes, int n, int x, int y, int val) {
    int num = get_val(boxes, n, x, y, 0);
    assert(num+1<MAX);
    set_val(boxes, n, x, y, num+1, val);
    set_val(boxes, n, x, y, 0, num+1);
}



//
//  benchmarking program
//
int main( int argc, char **argv )
{    
    if( find_option( argc, argv, "-h" ) >= 0 )
    {
        printf( "Options:\n" );
        printf( "-h to see this help\n" );
        printf( "-n <int> to set the number of particles\n" );
        printf( "-o <filename> to specify the output file name\n" );
        return 0;
    }

    int n = read_int( argc, argv, "-n", 1000 );

    char *savename = read_string( argc, argv, "-o", NULL );
    
    FILE *fsave = savename ? fopen( savename, "w" ) : NULL;
    particle_t *particles = (particle_t*) malloc( n * sizeof(particle_t) );
    set_size( n );
    init_particles( n, particles );

    // Initialize Boxes! woo!
    double size = 0.0;
    size = sqrt(density*n);
    int numBoxes =(int) ceil(size/cutoff);
    int *boxes = init_ds(numBoxes);


    //  simulate a number of time steps
    double simulation_time = read_timer( );
    for( int step = 0; step < NSTEPS; step++ )
    {

        // update boxes 
        for( int i = 0; i < n; i++){
            insert_ds(boxes, numBoxes, (int)(particles[i].x/cutoff),
                                (int)(particles[i].y/cutoff), i);
        }

        //  compute forces
        for( int i = 0; i < n; i++ )
        {
            particles[i].ax = particles[i].ay = 0;
            int currBox_x = (int) (particles[i].x/cutoff);
            int currBox_y = (int) (particles[i].y/cutoff);            

            for(int I = max(0,1-currBox_x); I < 3 && currBox_x-1+I<numBoxes; I++) {
                for(int J = max(0,1-currBox_y); J < 3 && currBox_y-1+J<numBoxes; J++){

                    int t = get_val(boxes, numBoxes, currBox_x-1+I, currBox_y-1+J, 0); 
                    for(int iter=1; iter <= t; iter++) {
                            apply_force( particles[i], particles[get_val(boxes, numBoxes, currBox_x-1+I, currBox_y-1+J, iter)]);
                    }

                }
            }
        }

        
        // clean boxes
        clean_ds(boxes, numBoxes); 

        //  move particles
        for( int i = 0; i < n; i++ ) 
            move( particles[i] );


        //  save if necessary
        if( fsave && (step%SAVEFREQ) == 0 )
            save( fsave, n, particles );
    }
    simulation_time = read_timer( ) - simulation_time;
    
    printf( "n = %d, simulation time = %g seconds\n", n, simulation_time );

    free( particles );
    if( fsave )
        fclose( fsave );
    
    return 0;
}
